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Abstract 

The well-known Stormer-Cowell class of linear k-step methods 
for the solution of second order initial value problems suffer 
from orbital instability. The solution of a problem describing a 
uniform motion in a circular orbit, spirals inward. Several 
modifications were suggested, some require the a-priori knowledge 
of the frequency . Here we develop a method based on the conser- 
vation of energy and momentum. This method overcomes the 
aforementioned difficulty. 
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1 . 



Introduction 



In this paper we develop a method for the numerical solution 
of the equations of motion of an object acted upon by several 
gravitational masses. In general, the motion can be described by 
a special class of second order differential equations, namely 



y"(x) = f(x,y(x)) . 



( 1 ) 



There exist methods of Runge-Kutta type which tackle this 
problem directly (Collatz [8, p. 61], de Vogelaere [23], Scraton 
[20] and others) . Also, linear multistep methods of the form 



k 

l a j¥n+j 

j = 0 



h2 



k* 

I 

j = 0 



& j f n+j 



(2) 



exist. See, for example, Henrici [11, p. 289], Lambert [13, p. 
252] and others. One of the authors has developed hybrid methods 
for such classes (Neta and Lee [15] , Neta [16, 18] ) . 

The direct application of methods of class (2) to problem 
(1), rather than the application of a conventional linear 
multistep method to an equivalent first-order system is usually 
recommended (Ash [1]). 

Special methods based on a-priori knowledge of the frequency 
were developed by Bettis [3], Steifel and Bettis [22], Gautschi 
[10], Neta and Ford [17], Neta [19], van der Houwen and Sommeijer 
[12], Lyche [14] and Sommeijer et al [21]. 
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To avoid the use of the frequency other methods were 
developed satisfying the so-called P-stability. This idea of P- 
stability was introduced by Lambert and Watson [13]. They showed 
that such methods are necessarily implicit and cannot have order 
greater than two. Chawla [5] and Cash [4] independently showed 
that this order-barrier can be crossed over by considering hybrid 
two-step methods. Other P-stable methods were developed by 
Costabile & Costabile [9], Chawla [6], Chawla and Rao [7] and 
others . 

Our idea here is to develop a new method to approximate the 
solution of the two body problem. This scheme will conserve both 
the energy per unit mass and the specific angular momentum of the 
system. 

In the next section the method is developed for the 
perturbation free flight. Numerical experiments with the scheme 
will be described in the last section. 
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2. Development of the Method 

In this section we describe our new scheme to solve the 
following system of equations 



d 2 x 

dt 2 



, 2 , 2, 3/2 
(x +y ) 



(3) 




, 2 2 > 3/2 

(x +y ) 



(4) 



where k = 3.9877 10 14 m 3 sec -2 is the gravitational parameter. 

It is well known that Cowell's method gives a numerical 
solution that spirals inward. Encke's method improves this 
result but requires more work [2]. 

It can be easily shown (see e.g., Bate et al. [2]) that the 
energy and momentum are conserved for a perturbation free flight. 
The conservation of these quantities will be used to obtain the 
fast changing variable. It is thus important to rewrite the 
system in polar coordinates 



d 2 r 




/ d S \ 2 

r( dt> 



k_ 

2 



(5) 



d 2 e _ 2 dr d9 (6) 

, , 2 r dt dt 

dt 



The radius r does not change much in time and thus we can use any 
method? here we use a second order Taylor series method. The 
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value of 9 corresponding to r will be obtained by integrating the 
relation (conservation of energy per unit mass) 



( d£) 2 + (r — ) 2 

v dt 1 dt ' 



<|f I ) 2 + (r(0) J|| ) 2 



t=0 



dt 



t=0 



(7) 



To be more specific, our method approximates r using 



r(t i+1 ) = r ( fc i) + r(ti) (At) + |r(t i )(At) 2 (8) 



where 



r(t i+ i) = r(t L ) + r ( i ) (At) (9) 

and r(t) is obtained from the differential equation. In a 
similar fashion we find 



0(t i + l> 



6 ( t i ) 



e ( t^) (At) 



(t ± ) (At) 



( 10 ) 



where 



6(t.) = \j (C + k/r(t i )) 2 - r 2 (t i )/r(t.) (11) 

and 9(t^) is obtained by differentiating the above relation. 

The results of this method were compared to Taylor series 
method for both r and 9 where now 9 is evaluated from an equation 
similar to (9) . 
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Remark: It should be clear that this idea can be used for any 
conservative system. One uses the conservation of some quantity 
to obtain an approximate value for the faster changing component. 
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3 . Numerical Experiments 

In this section we solve the equations of motion (3) -(4) 
combined with the initial values 



r ( 0) = R + 10 5 , 

r(0) = 0 , (12) 

e(0) = 0 , 

6(0) = 8000/r (0) , 

where R = 6.371*10 6 is the radius of earth. 

In the first table we compare the results of our method with 
Taylor series using t = 1 sec. As can be seen in the table, the 
value of r at the end of each period is nearly constant with our 
method but decreases as the number of periods increases with 
Taylor series method. 



TABLE 1 



At the 


7 




Energy per 


Specific angula 


end of 


r (10 


____ 


unit mass 


momentum 


Period 


Taylor 


New 


(10 8 m 2 /sec 2 ) 


(10-*- 1 m 2 /sec) 


0 


. 6471 


.6471 


-.2962 


. 5277 


5 


. 6320 


. 6476 


-.3019 


.5127 


10 


.6173 


.6482 


-.3077 


. 5078 


15 


.6030 


.6487 


-.3135 


. 5029 


20 


. 5890 


. 6493 


-.3194 


.4982 


30 


. 5622 


. 6503 


-.3313 


.4888 


40 


. 5367 


. 6513 


-.3433 


. 4798 


50 


. 5123 


. 6523 


-.3557 


.4710 


60 


. 4891 


.6533 


-.3683 


.4624 


70 


.4667 


.6543 


-.3811 


. 4540 
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The energy per unit mass E and the specific angular momentum are 
given by 



E = i(r 2 + (r0) 2 ) - | , 



P 



r 0 - r0 



(13) 



In our example both the energy and momentum are conserved. 

In the following three figures the radius r, the energy per 
unit mass E and the specific angular momentum P were plotted at 
the end of every 10 periods. It is clear that the relative 
change in r values after 70 periods using Taylor series method is 
more than 25 times larger compared to our new method. The energy 
and momentum were kept constant by the new scheme whereas the 
Taylor series method sho\*s a relative loss of 28% in energy and 
12% in momentum at the end of 70 periods. 

Remark : Note that our method is explicit whereas the P-stable 

methods suggested by Lambert and Watson [13] are implicit. 
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